function [MovAverage_reduced, MovAverage_structural] = zFC_RewriteVAR(Autoregressive, B, forhowlong)

% -----------------------------------------------------%
% DESCRIPTION OF THE FUNCTION
% This function takes a VAR written in auterregressive form and computes
% the moving average representations, both for the structural form and for
% the reduced form, given a B matrix that converts one into the other

% More precisely, the initial VAR model is written as
%
% y_t = Autoregressive y_t + B s_t
%     = A_1 y_t-1 + A_2 y_t-2 + ... + A_p y_t-p + B s_t
%
% where "Autoregressive" includes several lags, i.e. it is A(L), and "B"
% maps the structural shocks "s_t" into reduced form shocks. For
% "forhowlong" lags, the function computes the moving average reduced form
% representation
%
% y_t = MovAverage_reduced r_t
%     = r_t + C_1 r_t-1 + C_2 r_t-2 + ...
%
% where "MovAverage_reduced" includes forhowlong-1 lags (i.e. not all the
% infinite lags of the MA). Then, it computes the structural version of the
% moving average representation, i.e. 
%
% y_t = MovAverage_structural r_t
%     = D_0 s_t + D_1 s_t-1 + D_2 s_t-2 + ...
%
% where D_i = C_i * B
% -----------------------------------------------------%

% Companion form
k = size(Autoregressive,1);
p = size(Autoregressive,2)/k;
B_step = [Autoregressive; eye((p-1)*k), zeros((p-1)*k,k)];

% Moving average, reduced form representation
C_step = NaN*ones(k*p,k*p,forhowlong);
C_step(:,:,1) = eye(k*p);
C_step(:,:,2) = B_step;
for h = 3:forhowlong
   C_step(:,:,h) = C_step(:,:,h-1)*B_step;
end

MovAverage_reduced  = NaN*ones(k,k,forhowlong);
for h = 1:forhowlong
   MovAverage_reduced(:,:,h) = C_step(1:k,1:k,h);
end


%% MOVING AVERAGE STRUCTURAL FORM REPRESENTATION 
MovAverage_structural  = NaN*ones(k,k,forhowlong);
for h = 1:forhowlong
   MovAverage_structural(:,:,h) = MovAverage_reduced(:,:,h)*B;
end

end